Fritz: All Fritz/All Fritz.zip/All Fritz/FILES/PROGMISC/PCSSP.LZH/PC-SSP.ZIP/MATINSL.ZIP/DMFGR.FOR

Text File | 1985-11-29 | 5KB | 202 lines

C
C     ..................................................................
C
C        SUBROUTINE DMFGR
C
C        PURPOSE
C           FOR A GIVEN M BY N MATRIX THE FOLLOWING CALCULATIONS
C           ARE PERFORMED
C           (1) DETERMINE RANK AND LINEARLY INDEPENDENT ROWS AND
C               COLUMNS (BASIS).
C           (2) FACTORIZE A SUBMATRIX OF MAXIMAL RANK.
C           (3) EXPRESS NON-BASIC ROWS IN TERMS OF BASIC ONES.
C           (4) EXPRESS BASIC VARIABLES IN TERMS OF FREE ONES.
C
C        USAGE
C           CALL DMFGR(A,M,N,EPS,IRANK,IROW,ICOL)
C
C        DESCRIPTION OF PARAMETERS
C           A      - DOUBLE PRECISION GIVEN MATRIX WITH M ROWS
C                    AND N COLUMNS.
C                    ON RETURN A CONTAINS THE TRIANGULAR FACTORS
C                    OF A SUBMATRIX OF MAXIMAL RANK.
C           M      - NUMBER OF ROWS OF MATRIX A.
C           N      - NUMBER OF COLUMNS OF MATRIX A.
C           EPS    - SINGLE PRECISION TESTVALUE FOR ZERO AFFECTED BY
C                    ROUNDOFF NOISE.
C           IRANK  - RESULTANT RANK OF GIVEN MATRIX.
C           IROW   - INTEGER VECTOR OF DIMENSION M CONTAINING THE
C                    SUBSCRIPTS OF BASIC ROWS IN IROW(1),...,IROW(IRANK)
C           ICOL   - INTEGER VECTOR OF DIMENSION N CONTAINING THE
C                    SUBSCRIPTS OF BASIC COLUMNS IN ICOL(1) UP TO
C                    ICOL(IRANK).
C
C        REMARKS
C           THE LEFT HAND TRIANGULAR FACTOR IS NORMALIZED SUCH THAT
C           THE DIAGONAL CONTAINS ALL ONES THUS ALLOWING TO STORE ONLY
C           THE SUBDIAGONAL PART.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           GAUSSIAN ELIMINATION TECHNIQUE IS USED FOR CALCULATION
C           OF THE TRIANGULAR FACTORS OF A GIVEN MATRIX.
C           COMPLETE PIVOTING IS BUILT IN.
C           IN CASE OF A SINGULAR MATRIX ONLY THE TRIANGULAR FACTORS
C           OF A SUBMATRIX OF MAXIMAL RANK ARE RETAINED.
C           THE REMAINING PARTS OF THE RESULTANT MATRIX GIVE THE
C           DEPENDENCIES OF ROWS AND THE SOLUTION OF THE HOMOGENEOUS
C           MATRIX EQUATION A*X=0.
C
C     ..................................................................
C
      SUBROUTINE DMFGR(A,M,N,EPS,IRANK,IROW,ICOL)
C
C        DIMENSIONED DUMMY VARIABLES
      DIMENSION A(1),IROW(1),ICOL(1)
      DOUBLE PRECISION A,PIV,HOLD,SAVE
C
C        TEST OF SPECIFIED DIMENSIONS
      IF(M)2,2,1
    1 IF(N)2,2,4
    2 IRANK=-1
    3 RETURN
C        RETURN IN CASE OF FORMAL ERRORS
C
C
C        INITIALIZE COLUMN INDEX VECTOR
C        SEARCH FIRST PIVOT ELEMENT
    4 IRANK=0
      PIV=0.D0
      JJ=0
      DO 6 J=1,N
      ICOL(J)=J
      DO 6 I=1,M
      JJ=JJ+1
      HOLD=A(JJ)
      IF(DABS(PIV)-DABS(HOLD))5,6,6
    5 PIV=HOLD
      IR=I
      IC=J
    6 CONTINUE
C
C        INITIALIZE ROW INDEX VECTOR
      DO 7 I=1,M
    7 IROW(I)=I
C
C        SET UP INTERNAL TOLERANCE
      TOL=ABS(EPS*SNGL(PIV))
C
C        INITIALIZE ELIMINATION LOOP
      NM=N*M
      DO 19 NCOL=M,NM,M
C
C        TEST FOR FEASIBILITY OF PIVOT ELEMENT
    8 IF(ABS(SNGL(PIV))-TOL)20,20,9
C
C        UPDATE RANK
    9 IRANK=IRANK+1
C
C        INTERCHANGE ROWS IF NECESSARY
      JJ=IR-IRANK
      IF(JJ)12,12,10
   10 DO 11 J=IRANK,NM,M
      I=J+JJ
      SAVE=A(J)
      A(J)=A(I)
   11 A(I)=SAVE
C
C        UPDATE ROW INDEX VECTOR
      JJ=IROW(IR)
      IROW(IR)=IROW(IRANK)
      IROW(IRANK)=JJ
C
C        INTERCHANGE COLUMNS IF NECESSARY
   12 JJ=(IC-IRANK)*M
      IF(JJ)15,15,13
   13 KK=NCOL
      DO 14 J=1,M
      I=KK+JJ
      SAVE=A(KK)
      A(KK)=A(I)
      KK=KK-1
   14 A(I)=SAVE
C
C        UPDATE COLUMN INDEX VECTOR
      JJ=ICOL(IC)
      ICOL(IC)=ICOL(IRANK)
      ICOL(IRANK)=JJ
   15 KK=IRANK+1
      MM=IRANK-M
      LL=NCOL+MM
C
C        TEST FOR LAST ROW
      IF(MM)16,25,25
C
C        TRANSFORM CURRENT SUBMATRIX AND SEARCH NEXT PIVOT
   16 JJ=LL
      SAVE=PIV
      PIV=0.D0
      DO 19 J=KK,M
      JJ=JJ+1
      HOLD=A(JJ)/SAVE
      A(JJ)=HOLD
      L=J-IRANK
C
C        TEST FOR LAST COLUMN
      IF(IRANK-N)17,19,19
   17 II=JJ
      DO 19 I=KK,N
      II=II+M
      MM=II-L
      A(II)=A(II)-HOLD*A(MM)
      IF(DABS(A(II))-DABS(PIV))19,19,18
   18 PIV=A(II)
      IR=J
      IC=I
   19 CONTINUE
C
C        SET UP MATRIX EXPRESSING ROW DEPENDENCIES
   20 IF(IRANK-1)3,25,21
   21 IR=LL
      DO 24 J=2,IRANK
      II=J-1
      IR=IR-M
      JJ=LL
      DO 23 I=KK,M
      HOLD=0.D0
      JJ=JJ+1
      MM=JJ
      IC=IR
      DO 22 L=1,II
      HOLD=HOLD+A(MM)*A(IC)
      IC=IC-1
   22 MM=MM-M
   23 A(MM)=A(MM)-HOLD
   24 CONTINUE
C
C        TEST FOR COLUMN REGULARITY
   25 IF(N-IRANK)3,3,26
C
C        SET UP MATRIX EXPRESSING BASIC VARIABLES IN TERMS OF FREE
C        PARAMETERS (HOMOGENEOUS SOLUTION).
   26 IR=LL
      KK=LL+M
      DO 30 J=1,IRANK
      DO 29 I=KK,NM,M
      JJ=IR
      LL=I
      HOLD=0.D0
      II=J
   27 II=II-1
      IF(II)29,29,28
   28 HOLD=HOLD-A(JJ)*A(LL)
      JJ=JJ-M
      LL=LL-1
      GOTO 27
   29 A(LL)=(HOLD-A(LL))/A(JJ)
   30 IR=IR-1
      RETURN
      END